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ABSTRACT 

Asteroids residing in the first-order mean motion resonances with Jupiter hold impor- 
tant information about the processes that set the final architecture of giant planets. 
Here we revise current populations of objects in the J2/1 (Hecuba-gap group), J3/2 
(Hilda group) and J4/3 (Thule group) resonances. The number of multi-opposition as- 
teroids found is 274 for J2/1, 1197 for J3/2 and 3 for J4/3. By discovering a second and 
third object in the J4/3 resonance, (186024) 2001 QG207 and (185290) 2006 UB219, 
this population becomes a real group rather than a single object. Using both hierarchi- 
cal clustering technique and colour identification we characterise a coUisionally-born 
asteroid family around the largest object (1911) Schubart in the J3/2 resonance. There 
is also a looser cluster around the largest asteroid (153) Hilda. Using iV-body numeri- 
cal simulations we prove that the Yarkovsky effect (infrared thermal emission from the 
surface of asteroids) causes a systematic drift in eccentricity for resonant asteroids, 
while their semimajor axis is almost fixed due to the strong coupling with Jupiter. 
This is a different mechanism from main belt families, where the Yarkovsky drift af- 
fects basically the semimajor axis. We use the eccentricity evolution to determine the 
following ages: (1.7 ± 0.7) Gyr for the Schubart family and > 4 Gyr for the Hilda fam- 
ily. We also find that collisionally-born clusters in the J2/1 resonance would efficiently 
dynamically disperse. The steep size distribution of the stable population inside this 
resonance could thus make sense if most of these bodies are fragments from an event 
older than ~ 1 Gyr. Finally, we test stability of resonant populations during Jupiter's 
and Saturn's crossing of their mutual mean motion resonances. In particular we find 
primordial objects in the J3/2 resonance were efficiently removed from their orbits 
when Jupiter and Saturn crossed their 1:2 mean motion resonance. 
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1 INTRODUCTION 

Populations of asteroids in the Jovian first order mean mo- 
tion resonances -J2/1, J3/2 and J4/3- are closely linked to 
the orbital evolution of the giant planets. This is because of 
their orbital proximity to Jupiterlj Stability or instability of 
these asteroid populations directly derives from the orbital 
configuration of the giant planets. As such it is also sensitive 
on the nature and amount of Jupiter's migration and other 
finer details of its dynamics. As a result, the currently ob- 
served asteroids in the Jovian first order resonances contain 
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residing in the J3/2 and J4/3 resonances, immediately attracted 
attention of astronomers by vastly extending asteroid zone toward 
giant planets and by their ability to apparently approach Jupiter 
near aphelia of their orbits (e.g., Kiihnert 1876; Krueger 1889). 



valuable information about the early evolution of planets 
and, if correctly understood and properly modelled, they 
may help to constrain it. 



Apart from the Trojan clouds (not studied in this paper) 
the largest known population in the Jovian mean motion 
resonances occupies the J3/2 resonance, and is frequently 
called the Hilda group. It was carefully studied in a par- 
allel series of papers by Schubart and Dahlgren and col- 
laborators during the past few decades. Schubart (1982a,b, 
1991) analysed short-term dynamics of Hilda-type orbits 
and introduced quasi-constant orbital parameters that al- 
lowed their first classification. While pioneering, Schubart 's 
work had the disadvantage of having much smaller sample of 
known asteroids and computer power than today. Dahlgren 
& Lagerkvist (1995) and Dahlgren et al. (1997, 1998, 1999) 
conducted the first systematic spectroscopic and rotation- 
rate investigation of Hildas. They found about equal abun- 
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dance of D- and P-type asteroid^ and suggested spectral- 
size correlation such that P-types dominate large Hildas and 
D-type dominate smaller Hildas. They also suggested small 
Hildas have large lightcurve amplitudes, as an indication of 
elongated or irregular shape, and that the distribution of 
their rotation rates is non-Maxwellian. Further analysis us- 
ing the Sloan Digital Sky Survey (SDSS) data however does 
not support significant dominance of either of the two spec- 
tral types for small sizes and indicates about equal mix of 
them (Gil-Hutton & Brunini 2008; see also below). Smaller 
populations of asteroids in the J2/1 and J4/3 received com- 
paratively less observational effort. 

Since the late 1990s powerful-enough computers allowed 
a more systematic analysis of fine details of the longer-term 
dynamics in the Jovian first order resonances. Ferraz-Mello 
& Michtchenko (1996) and Ferraz-MeUo et al. (1998a,b) de- 
termined that asteroids in the J2/1 resonance can be very 
long-lived, possibly primordial, yet their motion is compar- 
atively more chaotic than those in the J3/2 resonance. The 
latter paper showed that commensurability between the li- 
bration period and the period of Jupiter's and Saturn's 
Great Inequality might have played an important role in 
depletion of the J2/1 resonance. This would have occurred 
when both giant planets were farther from their mutual 2:5 
mean motion configuration in the past. A still more complete 
analysis was obtained by Nesvorny & Ferraz-Mello (1997) 
who also pointed out that the J4/3 resonance stable zone is 
surprisingly void of asteroids, containing only (279) Thule. 
Roig et al. (2001) and Broz et al. (2005) recently revised the 
population of asteroids in the J2/1 resonance and classified 
them into several groups according to their long-term orbital 
stability. While the origin of the unstable resonant popula- 
tion was successfully interpreted using a model of a steady- 
state flow of main belt objects driven by the Yarkovsky semi- 
major axis drift, the origin of the long-lived asteroids in the 
J2/1 remains elusive. Population of Hildas and Thule was 
assumed primordial or captured by an adiabatic migration 
of Jupiter (e.g.. Franklin et al. 2004). 

It has been known for some time that the current con- 
figuration of giant planets does not correspond to that at 
their birth. However, a new momentum to that hypothesis 
was given by the so called Nice model (Tsiganis et al. 2005; 
Morbidelh et al. 2005; Gomes et al. 2005). The Nice model 
postulates the initial configuration of the giant planets was 
such that Jupiter and Saturn were interior of their mutual 
1:2 mean motion resonance (see also Morbidelli et al. 2007). 
The event of crossing this resonance had a major influence 
on the final architecture of giant planets and strongly influ- 
enced structure of small-bodies populations in the Solar sys- 
tem. Morbidelli et al. (2005) showed that the population of 
Jupiters Trojan asteroids was destabilised and re-populated 
during this phase. In what follows we show that, within the 
Nice model, the same most probably occurs for populations 
of asteroids in the J3/2 and J4/3 resonances. 

The paper is organised as follows: in Section [2] we re- 
vise information about the current populations of asteroids 
in the Jovian first order resonances. We use an up-to-date 
AstOrb database of asteroid orbits from the Lowell Observa- 
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tory (ftp.lowell.edu) as of September 2007 and eliminate 
only single-opposition cases to assure accurate orbital infor- 
mation. 

In SectionOwe apply clustering techniques and extract 
two families of asteroids on similar orbits in the J3/2 reso- 
nance. We strengthen their case with an additional colour 
analysis using the SDSS broadband data. We model the 
long-term orbital evolution of these families and estimate 
their ages on the basis of Yarkovsky-driven dispersion in ec- 
centricity. 

In Section 13] we determine an orbital stability of the pu- 
tative primordial populations of planetesimals in the Jovian 
first-order resonances. We show that those in the J3/2 and 
J4/3 are very efficiently eliminated when Jupiter and Sat- 
urn cross their mutual 1:2 mean motion resonance. We also 
determine the removal rate of very small resonant asteroids 
due to the Yarkovsky/YORP effects. 



2 CURRENT ASTEROID POPULATIONS IN 
THE JOVIAN FIRST ORDER RESONANCES 

Dynamics of asteroid motion in the Jovian first order reso- 
nances has been extensively studied by both analytical and 
numerical methods in the past few decades (e.g., Murray 
1986; Sessin & Bressane 1988; Ferraz-MeUo 1988; Lemaitre 
& Henrard 1990; Morbidelh & Moons 1993; Moons et al. 
1998; Nesvorny & Ferraz-MeUo 1997; Roig et al. 2002; 
Schubart 2007). In what follows we review a minimum in- 
formation needed to understand our paper, referring an in- 
terested reader to the literature mentioned above for more 
insights. 

In the simplest framework of a circular restricted planar 
three-body problem (Sun- Jupiter-asteroid) the fundamental 
effects of the resonant dynamics is reduced to a one-degree 
of freedom problem defined by a pair of variables (E, a). For 
J(p -f- l)/p resonance (p = 1, 2 and 3 in our cases) we have 

E = (l-\/l-e2) , (1) 

a = {p+l)X' -pX-zu , (2) 

where a is the semimajor axis, e the eccentricity, zu the lon- 
gitude of pericentre and A the mean longitude in orbit of the 
asteroid, and A' is the mean longitude in orbit of Jupiter. 

If the asteroid motion is not confined into the orbital 
plane of the planet, we have an additional pair of resonant 
variables (Ez,^^) such that 

E, = 2Va(l-e2)sin2 ^ , (3) 

a, = (p+l)A'-pA-!^, (4) 

where i denotes the inclination of asteroids orbit and SI the 
longitude of its node. Remaining still with the simple av- 
eraged model, orbital effects with shorter periods are ne- 
glected, the motion obeys an integral of motion A'' given by 

iV = (^^^- \/l-e2cosi^ . (5) 

Because of this integral of motion, variations of E imply 
oscillations of both a and e. 

The two-degree of freedom character of the resonant 
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motion prevents integrability. However, as an approxima- 
tion we may introduce a hierarchy by noting that pertur- 
bation described by the (E, a) variables is larger than that 
described by the (E^jCTz) terms (e.g., Moons et al. 1998). 
This is usually true for real resonant asteroids of interest. 
Only the angle a librates and cr^ circulates with a very long 
period. The (Ez,(Jz) dynamics thus produces a long-period 
perturbation of the (E, cr) motion. 

Within this model the minimum value of E in one reso- 
nant cycle (typically several hundreds of years) implies a is 
minimum and e is maximum. These values do not conserve 
exactly from one cycle to another because the (Ez,az) mo- 
tion produces small oscillations. Since T. + 'E^ — N = —^/H/p 
one needs to wait until E^ reaches maximum over its cycle 
to attain 'real' minimum of a values and 'real' maximum 
of e values over a longer time interval. From (|3]) we note 
the maximum of Ez occurs for the maximum of i variations. 
This situations occurs typically once in a few thousands of 
years. In an ideal situation, these extremal values of (a, e, i) 
would be constant and may serve as a set of proper orbital 
elements. 

The motion of real asteroids in the Solar system is fur- 
ther complicated by Jupiter having non-zero and oscillat- 
ing value of eccentricity. This brings further perturbations 
(e.g., Ferraz-Mello 1988; Sessin & Bressane 1988 for a sim- 
ple analytic description) and sources of instability inside 
the resonance. Despite the non-integrability, we follow Roig 
et al. (2001) and introduce pseudo-proper orbital elements 
(op, Bp, sin ip) as the osculating elements (a,e,sini) at the 
moment, when the orbit satisfies the condition: 

a = OA^<OAzu~zu' = OAn-n' = 0, (6) 
dt . V / 

where uj' and Q' denote the longitude of pericentre and the 
longitude of node of Jupiter. As above, when ((6]) holds the 
osculating orbital elements are such that a attains minimum, 
e attains maximum and i attains maximum. Numerical ex- 
periments show that with a complete perturbation model 
and a finite time-step it is difficult to satisfy all conditions 
of ((B]) simultaneously. Following Roig et al. (2001) we thus 
relax ((6]) to a more practical condition 

|cr! < 5° A ^ < OA Icc-cc'l < 5° . (7) 

Because this condition is only approximate, we numerically 
integrate orbits of resonant asteroids for 1 Myr, over which 
the pseudo-proper orbital elements are recorded. We then 
compute their mean value and standard deviation, which is 
an expression of the orbital stability over that interval of 
time. 

In the case of the J3/2 and J4/3 resonances, we use 
the condition Q with a different sign ^ > and, more- 
over, we apply a digital filter (denoted as A in Quinn et al. 
(1991), using 1 yr sampling and a decimation factor of 10) to 
a{t). This intermediate stage serves to suppress oscillations 
faster than the libration period. The different sign of ^ 
just means our pseudo-proper orbital elements correspond 
to maximum value of a and minimum values of e and i, in 
order to allow more direct comparison with previous analy- 
ses. 

Aside to this short-term integration we perform long- 
term runs to determine the stability of a particular reso- 
nant orbit. With this aim we conduct integrations spanning 



4Gyr for all resonant asteroids. Because of the inherent un- 
certainty in the initial conditions (orbital elements at the 
current epoch), we perform such integration for the nomi- 
nal orbit and 10 clones that randomly span the uncertainty 
ellipsoid. We then define dynamical lifetime of the orbit as 
the median of time intervals, for which the individual clones 
stayed in the resonance. 

All integrations are performed using the Swift pack- 
age (Levison & Duncan 1994), slightly modified to include 
necessary on-line digital filters and a second order symplec- 
tic integrator (Laskar & Robutel 2001). Most of numeri- 
cal simulations take into account gravitational interactions 
only, but in specific cases - and when explicitly mentioned 
- we include also Yarkovsky (thermal) accelerations. In this 
case we use an implementation described in detail by Broz 
(2006). Our simulations include 4 outer planets. We mod- 
ify the initial conditions of the planets and asteroids by a 
barycentric correction to partially account for the infiuence 
of the terrestrial planets. The absence of the terrestrial plan- 
ets as perturbers is a reasonable approximation in the outer 
part of the main belt and for orbits with e < 0.8 in gen- 
eral. We nevertheless checked the short-term computations 
(determination of pseudo-proper resonant elements) using a 
complete planetary model and noticed no significant differ- 
ence in results. The second order symplectic scheme allows 
us to use a time step of 91 days. 

2.1 Hecuba-gap group 

In order to determine, which objects are located in the J2/1 
mean motion resonance, we first extracted orbits from the 
AstOrb database with osculating orbital elements in a broad 
box around this resonance (see, e.g., Roig et al. 2001 for a 
similar procedure). We obtained 7139 orbits, which we nu- 
merically integrated for 10 kyr. We recorded and analysed 
behaviour of the resonance angle a — 2X' — X ~ -ou from 
Eq. ((2}. Pericentric librators, for which a oscillates about 
0°, were searched. We found 274 such cases; this extends 
the previous catalogue of Broz et al. (2005) almost twice. 
The newly identified resonant objects are mainly asteroids 
discovered or recovered after 2005 with accurate enough or- 
bits. We disregard from our analysis asteroids at the border 
of the resonance, for which a{t) exhibits alternating peri- 
ods of libration and circulation, and also those asteroids for 
which a oscillates but are not resonant anyway (A^ ^ 0.8 
in Eq. (O; see, e.g., Morbidelli & Moons 1993). The latter 
reside on low-eccentricity orbits in the main asteroid belt 
adjacent to the J2/1 resonance. 

We conducted short- and long-term integrations of the 
resonant asteroids as described above. They allowed us to 
divide the population into 182 long-lived asteroids (with the 
median dynamical lifetime longer than 70 Myr, as defined 
in Broz et al. 2005) and 92 short-lived asteroids (the life- 
time shorter than 70 Myr) , see Fig. Among the short- 
lived objects we found 14 have dynamical lifetimes even less 

Our results for both J2/1 and J3/2 resonances 
are summarised in tables available through a website 
http://sirrah.troja.mff.cmii.cz/yarko-site/ (those for 
J4/3 bodies are given in Table [TJ. These contain listing of all 
resonant asteroids, their pseudo-proper orbital elements with 
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Figure 1. A distribution of the median dynamical lifetimes for 
objects in the J2/1 resonance. A division to several groups is 
denoted: extremely unstable objects (tj2/i 2 Myr), short-lived 
objects (tj2/i ^ 70 Myr) and long-Uved objects (Griquas and 
Zhongguos, tj2/i > 70 Myr). 



than 2 Myr and we call them extremely unstable. Broz et al. 
(2005) suggested the unstable orbits in the J2/1 resonance 
are resupplied from the adjacent main belt due to a per- 
manent flux driven by the Yarkovsky force, the extremely 
unstable objects are most probably temporarily captured 
Jupiter-family comets. The origin of the long-lived popula- 
tion in this resonance is still not known. 

Figure [2] shows the pseudo-proper orbital elements of 
the J2/1 asteroids projected onto the (a^, e^) and (op, sinip) 
planes. Our data confirm that the unstable population of 
J2/1 asteroids populates the resonance outskirts near its 
separatrix, where several secular resonances overlap and 
trigger chaotic motion (e.g., Morbidelli & Moons 1993; 
Nesvorny & Ferraz-Mello 1997; Moons et al. 1998). At low- 
eccentricities the chaos is also caused by an overlap on 
numerous secondary resonances (e.g., Lemaitre & Henrard 
1990). Two 'islands' of stability -A and B- harbour the long- 
lived population of bodies. The high-inclination island A, 
separated from the low-inclination island B by the ui^ sec- 
ular resonance, is much less populated. Our current search 
identifies 9 asteroids in the island A. The origin of the asym- 
metry in A/B islands is not known, but since the work of 
Michtchenko & Ferraz-Mello (1997) and Ferraz-MeUo et al. 
( 1998a, b) it is suspected to be caused by instability due to 
the libration period commensurability with the forcing terms 
produced by the Great Inequality. 

The size- frequency distribution of objects of a popula- 
tion is an important property, complementing that of the 
orbital distribution. Figure |3] shows cumulative distribution 
N{<H) of the absolute magnitudes H for bodies in the 
J2/1 (and other Jovian first order resonances as well). In 
between H — 12 mag and 14.5 mag (an approximate com- 
pleteness limit; R. Jedicke, personal communication) it can 
be matched by a simple power- law N{<H) oc 10^-^, with 
7 = (-1-0.70 ± 0.02)0 We thus confirm that the J2/1 pop- 



standard deviations, their dynamical residence time and some 
additional information. 

^ This is equivalent to a cumulative size distribution law 
N{>D) oc with a = —5-f = (—3.5 it 0.1), assuming all bodies 
have the same albedo. 
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Figure 2. Pseudo-proper orbital elements for the 247 objects 
in the J2/1 resonance projected onto the planes of semimajor 
axis ap vs eccentricity Cp (top) and semimajor axis ap vs sine of 
inclination sinip (bottom). Bars are standard deviations of the 
elements derived from 1 Myr numerical integration. Position of 
several secular resonances embedded in J2/ 1 is shown in the upper 
panel. The unstable population of asteroids (crosses) occupies 
the region of their overlap; the stable population (full circles) 
occupies two distinct zones —A and B— of low-eccentricity and 
low-inclination orbits (e.g., Nesvorny & Ferraz-Mello 1997). The 
population of marginally stable asteroids (open squares) resides 
in region adjacent to the unstable borders of the resonance or 
near the bridge over the stable regions associated with the i/ig 
secular resonance. 



ulation is steeper than it would correspond to a standard 
coUisionally evolved system (e.g., Dohnanyi 1969; O'Brien 
& Greenberg 2003) with 7 — +0.5. The same result holds 
for both the short- and long-lived sub-populations in this 
resonance separately. 

Albedos of J2/1 bodies are not known, except for 
(1362) Griqua for which Tedesco et al. (2002) give pv = 
0.067. The surroundings main-belt population has an aver- 
age Pv ~ 0.05. For sake of simplicity we convert absolute 
magnitudes to sizes using this averaged value when needed. 
For instance in Fig. |4]we show a zoom on the long-lived pop- 
ulation of objects in the J2/1 resonance with symbol size 
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Figure 3. Cumulative distributions N{<H) of the absolute mag- 
nitude H values for population of asteroids in the Jovian first 
order mean motion resonances: (i) J2/1 (dashed curve), (ii) J3/2 
(solid curve), and (iii) J4/3 (dotted curve). The straight lines 
show best-fit approximations N{<H) oc IC'^ with the values of 
7 indicated by the corresponding label. The fit matches N{<H) 
for H in the interval (12, 15) mag for J2/1 and (10.5, 14.5) mag for 
J3/2; no such approximation is available for J4/3 where only three 
objects arc currently known. The H values where the straight line 
approximations level off from the data roughly correspond to the 
completeness limit of the population (R. Jedicke, personal com- 
munication). For sake of a rough comparison, the upper abscissa 
gives an estimate of sizes for the albedo value py = 0.05, average 
of the outer belt population. 



weighted by the estimated size of the body. We note large 
objects are located far from each other and they are quite 
isolated — no small asteroids are in close surroundings. Both 
these observations suggest that the long-lived J2/1 popula- 
tion does not contain recently-born coUisional clusters. 
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Figure 4. A zoom on the (op, ep) and (op, sin ip) plots from Fig. [2] 
with relative size of the resonant asteroids indicated by size of the 
crosses. Note the large bodies, some of which are labelled, reside 
far from each other. 



2.2 Hilda group 

Because asteroids in the J3/2 constitute a rather isolated 
group, it is easy to select their candidates: we simply ex- 
tracted from the AstOrb database those asteroids with semi- 
major axis in between 3.8 AU and 4.1 AU. With that we 
obtained 1267 multi-opposition objects. We numerically in- 
tegrated these orbits for lOkyr and analysed the behaviour 
of the resonance angle a — 3A' — 2X — w. We obtained 1197 
cases for which a librates about 0° and which have A'' 0.44, 
a threshold of the resonance zone (e.g., Morbidelli & Moons 
1993); see Fig. [5] 

The long-term evolution of Hildas indicates that not 
all of them are stable over 4 Gyr, but 20 % escape earlier. 
A brief inspection of Fig. [6] shows, that the escapees are 
essentially asteroids located closer to the outer separatrix 
and exhibiting large amplitudes of librations. If the Hilda 
group has been constituted during the planetary formation 
some 4 Gyr ago, some non-conservative process must have 
placed these objects onto their currently unstable orbits. We 
suspect mutual collisions or gravitational scattering on the 



largest Hilda members might be the corresponding diffusive 
mechanisms. Small enough members might be also suscep- 
tible to the resonant Yarkovsky effect (see Sec. 14.21 and the 
Appendix A). 

Data in Fig. |3] confirm earlier findings that the Hilda 
group is characterized by an anomalously shallow size distri- 
bution. In between absolute magnitudes H = 10.5 mag and 
14.5 mag the cumulative distribution can be well matched 
by N{<H) oc 10^-^ with 7 = {+0.35 ± 0.02) only. 

Sub-populations among Hilda asteroids, namely two 
coUisional families, are studied in Section [S] 

2.3 Thule group 

In spite of a frequent terminology "Thule group" , asteroids 
in the J4/3 resonance consisted of a single object (279) Thule 
up to now. Nesvorny & Ferraz-Mello (1997) considered this 
situation anomalous because the extent of the stable zone of 
this resonance is not much smaller than that of the J3/2 res- 
onance (see also Franklin et al. 2004). In the same way, our 
knowledge about the low-e and low-i Thule-type stable or- 
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Figure 5. Pseudo-proper orbital elements for 1197 Hildas pro- 
jected onto tlie planes of semimajor axis ap vs eccentricity Sp 
(top) and semimajor axis ap vs sine of inclination sinip (bottom). 
Larger size of the symbol indicates larger physical size of the as- 
teroid. Because of Hildas orbital stability, the uncertainty in the 
pseudo-proper element values is typically smaller than the sym- 
bol size. Note a tight cluster around the proper inclination value 
sinip ~ 0.0505, lead by the largest asteroid (1911) Schubart, and 
a somewhat looser cluster around the proper inclination value 
sinip ~ 0.151, lead by the largest asteroid (153) Hilda. Both are 
discussed in more detail in Section [S] Solid line denotes the libra- 
tion centre of the J3/2 resonance. 



bits (e.g., a = 4.27 AU, e = 0.1 and i = 5°) should be obser- 
vationally complete at about magnitudes H = 12.5-13 mag 
(R. Jedicke, personal communication). A rough estimate also 
shows that even one magnitude in H beyond this complete- 
ness limit the Thule population should be known at ~ 10% 
completeness, leaving only about 90 % undiscovered popu- 
lation. We thus conclude that the objects in the magnitude 
range H = 9-13 mag are very likely missing in this reso- 
nance. Where does the existing population of small Thule- 
type asteroids begin? 

Our initial search in the broad box around the J4/3 
resonance detected only 13 objects. Six of them, includ- 
ing the well-known extinct comet (3552) Don Quixote (e.g., 
Weissman et al. 2002), are on typical orbits of Jupiter- family 
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Figure 6. Pseudo-proper semimajor axis vs eccentricity plot for 
asteroids in the J3/2 resonance. Bold crosses denote orbits, which 
escaped during the 4Gyr of evolution. Solid line is the position 
of the libration centre (the outer separatrix is located further to 
the right). 
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Figure 7. The histogram shows distribution of values of the ge- 
ometric albedo py determined by Tedesco et al. (2002) for aster- 
oids located inside the J3/2 resonance. Three individual values - 
(153) Hilda, (1911) Schubart and (279) Thule (a J4/3 asteroid) - 
are indicated on top. 



comets that happen to reside near this resonance with very 
high eccentricity and moderately high inclination. Two more 
are single-opposition objects and one has only poorly con- 
strained orbit, leaving us with (279) Thule and three ad- 
ditional candidate objects: (52007) 2002 EQ47, (186024) 
2001 QG207 and (185290) 2006 UB219. 

Figure |8] (top panels) shows short-term tracks of (279) 
Thule, (186024) 2001 QG207 and (185290) 2006 UB219 in res- 
onant variables \/2E (cos a, sin a) ~ (e cos cr, e sin cr) of the 
J4/3 resonance {a = 4A' — 3A — in this case). In all cases 
their orbits librate about the pericentric branch (a = 0°) of 
this resonance, although this is complicated -mainly in the 
low-eccentricity case of (279) Thule- by the forced terms due 
to Jupiter's eccentricity (see, e.g., Ferraz-Mello 1988; Sessin 
& Bressane 1988; Tsuchida 1990). The leftmost panel re- 
covers the 40°-50° libration of (279) Thule, determined pre- 
viously by Tsuchida (1990; Fig. 3). The other two smaller 
asteroids show librations with comparable amplitudes. The 
last object, (52007) 2002 EQ47, appears to reside on an un- 
stable orbit outside the J4/3 resonance. Our search thus 
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Table 1. Data on presently known population of asteroids residing in the J4/3 Jovian mean motion reso- 
nance (Thule group). Pseudo-proper orbital elements (ap, ep, sin jp) are given together with their standard 
deviations {Sap, Sep, S sinip) determined from a 1 Myr numerical integration. Up^rnax is the maximum libra- 
tion amplitude in the Sessin's (K,H) variables (see Fig. (Sjl, H is the absolute magnitude from the AstOrb 
catalogue and D is the estimated size using py = 0.04 geometric albedo (Tedesco et al. 2002). 



No. 


Name 


ap 


ep 


sinip 


Sap 


Sep 


S sin ip 


'^p, max 


H 


D 






[AU] 






[AU] 






[deg] 


[mag] 


[km] 
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4.2855 


0.119 


0.024 
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0.012 


0.003 


~50 


8.57 


126.6 


186024 


2001 QG207 


4.2965 


0.244 


0.042 
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0.014 


0.003 


25 


14.36 
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Figure 8. Top panels: orbits of (279) Thule [left], (186024) 2001 QG207 [middle] and (185290) 2006 UB219 [right] in resonant variables 
(e cos (J, e sin cr) of the J4/3 resonance. Here e is the eccentricity and a = AX' — 3A — ro, with A and A' the mean longitude in orbit of the 
asteroid and Jupiter and ro is the longitude of asteroid's pericentre (in all cases the osculating orbital elements are used). Each panel 
shows results of a short-term, 10 kyr numerical integration. In all three cases the orbits librate about the pericentric branch (tr = 0°) of 
the resonance. Perturbations due to Jupiter's eccentricity and its variations make the regular libration move in an epicyclic manner (e.g., 
Ferraz-Mello 1988; Sessin & Brassane 1988). Bottom panels: filtered resonant variables (e cos cr, esin cr), with short-period oscillations 
forced by Jupiter removed by digital filtering. They are similar to Sessin's {K, H) coordinates in Sessin & Brassane (1988) or Tsuchida 
(1990). In each of the cases the pericentric libration is clearly revealed. The maximum libration amplitude in these coordinates is denoted 
by o-p^rnax in Table [1] 



lead to the detection of two new asteroids in this resonance, 
increasing its population by a factor of threelfl 

Results of a long-term numerical integration of the nom- 
inal orbits plus 10 close clones, placed within an orbital un- 
certainty, reveal that the orbit of (279) Thule is stable over 
4Gyr, but the orbits of (186024) 2001 QG207 and (185290) 
2006 UB219 are partially unstable. They are not 'short- 
lived' but 45 % and 60 % of clones, respectively, escaped be- 
fore 4Gyr. Figure [§] shows pseudo-proper semimajor axis vs 
time for nominal orbits and their clones of all J4/3 objects; 
the escaping orbits leave the figure before the simulation 



^ While our initial search used AstOrb catalogue from Septem- 
ber 2007, we repeated it using the catalogue as of June 2008. No 



additional J4/3 objects were found. 



was ended at 4 Gyr. We suspect similar non-conservative ef- 
fects as mentioned above for the 20% fraction of long-term- 
unstable Hildas to bring these two Thule members onto their 
marginally stable orbits. 

We would like to point out that it took more than a cen- 
tury from the discovery of (279) Thule (Palisa 1888; Krueger 
1889) until further objects in this resonance were finally dis- 
covered. This is because there is an anomalously large gap 
in size of these bodies: (279) Thule is 127 km in size with 
pv = 0.04 AU (e.g., Tedesco et al. 2002), while the estimated 
sizes of (186024) 2001 QG207 and (185290) 2006 UB219 for 
the same value of albedo are 8.9 km and 11.8 km only. It 
will be interesting to learn as much as possible about the 
Thule population in the H = 13-15 mag absolute magnitude 
range using future-generation survey projects such as Pan- 
STARRS (e.g., Jedicke et al. 2007). Such a completed pop- 
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Figure 9. Pseudo-proper semimajor axis ap vs time t plot for 
the 3 orbits of asteroids located inside tlie J4/3 resonance and 
10 close clones orbits (placed randomly within the uncertainty 
ellipsoids) for each of them. Thin lines denote the stable orbits 
and thick lines unstable, which escaped from the J4/3 before the 
completion of the simulation at 4Gyr. 



ulation may present an interesting constraint on the plan- 
etesimal size distribution 4Gyr ago. 
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Figure 10. The dependence of the minimum number of asteroids 
^mini to be considered a statistically significant cluster, on the 
cutoff velocity iicutoff (thick curve). The average number A^o and 
the maximum number max{Ni) of asteroids, which are closer than 
'f'cutoffi is shown by dashed and thin curves. All quantities are 
valid for the J3/2 population. The fact that max(Afi) is much 
larger than A^o and A'^min indicates a presence of a significant 
cluster (or clusters) among the Hilda group. 



3 COLLISIONAL FAMILIES AMONG HILDA 
ASTEROIDS 

Collisions and subsequent fragmentations are ubiquitous 
processes since planets formed in the Solar system. Because 
the characteristic dispersal velocities of the ejecta (as a rule 
of thumb equal to the escape velocity of the parent body) 
are usually smaller than the orbital velocity, the resulting 
fragments initially reside on nearby orbits. If the orbital 
chaoticity is not prohibitively large in the formation zone, 
we can recognise the outcome of such past fragmentations as 
distinct clusters in the space of sufficiently stable orbital ele- 
ments. More than 30 collisional families are known and stud- 
ied in the main asteroid belt (e.g., Zappala et al. 2002) with 
important additions in the recent years (e.g., Nesvorny et al. 
2002; Nesvorny, Vokrouhlicky & Bottke 2006). Similarly, col- 
lisional families have been found among the Trojan clouds of 
Jupiter (e.g., Milani 1993; Beauge & Roig 2001; Roig et al. 
2008), irregular satellites of Jupiter (e.g., Nesvorny et al. 
2003, 2004) and even trans-Neptunian objects (e.g.. Brown 
et al. 2007). Mean motion resonances, other than the Trojan 
librators of Jupiter, are typically too chaotic to hold sta- 
ble asteroid populations, or the populations were too small 
to enable search for families. The only remaining candidate 
populations are those in the Jovian first order resonances, 
with Hilda asteroids the most promising group. However, 
low expectations for an existence of collisional families likely 
de-motivated systematic search. Note, the estimated intrin- 
sic collisional probability of Hilda asteroids is about a factor 
3 smaller than in the main asteroid belt (e.g., Dahlgren 1998; 
DeirOro et al. 2001) and the population is more than two 
orders of magnitude smaller. 

In spite of the situation outlined above, Schubart 
(1982a, 1991) repeatedly noticed groups of Hilda-type as- 
teroids with very similar proper elements. For instance, in 
his 1991 paper he lists 5 members of what we call Schubart 
family below and pointed out their nearly identical values of 



the proper inclination. Already in his 1982 paper Schubart 
mentions a similarity of such clusters to Hirayama families, 
but later never got back to the topic to investigate this prob- 
lem with sufficient amount of data provided by the growing 
knowledge about the J3/2 populationlf] Even a zero order 
inspection of Fig. [S] in particular the bottom panel, implies 
the existence of two large clusters among the J3/2 popu- 
lation. In what follows we pay a closer analysis to both of 
them. 

We adopt an approach similar to the hierarchical- 
clustering method (HCM) frequently used for identification 
of the asteroid families in the main belt (e.g., Zappala et al. 
1990, 1994, 2002). In the first step of our analysis, we com- 
pute the number of bodies A^min which is assumed to con- 
stitute a statistically significant cluster for a given value 
of the cutoff velocity iicutoff - We use a similar approach to 
that of Beauge & Roig (2001): for all asteroids in the J3/2 
resonance we determine the number Ni{vcuto{!) of asteroids 
which are closer than Wcutoff- Then we compute the average 
value A'^o = Ni. According to Zappala et al. (1994), a cluster 
may be considered significant if A'^ > A^'min = A^o + 2\/ No- 
The plots A'o(wcutoff ) and the corresponding Afmin(t'cutoff ) for 
Hilda population are shown in Fig. 1101 We use a standard 
metric (di defined by Zappala et al., 1994), namely: 

5v = nap ^^(^) +2(<5ej,)2 + 2((5sinv)2 , (8) 

where (a^, Bp, sin ip) are 10 Myr averaged values of the res- 
onant pseudo-proper elements (we checked that our results 
practically do not depend on the width of this averaging 
interval) . 

^ Schubart lists 11 additional asteroids in the group on his web- 
site http : //www . rzuser . uni-heidelberg . de/~ s24/hilda . htm 
but again he does not go into details of their putative collisional 
origin. 
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Figure 11. A stalactite diagram computed for ttie J3/2 popula- 
tion (Hildas). Two prominent groupings, the Schubart family and 
the Hilda family, are indicated. Every group plotted here has at 
least 5 or A'^niin members, whichever is larger (see Fig. IIOII . 



Next, we construct a stalactite diagram for Hildas in 
a traditional way (e.g., Zappala et al. 1990): we start with 
(153) Hilda as the first central body and we find all bodies 
associated with it at Wcutofi = 300 m/s, using a hierarchical 
clustering method (HCM; Zappala et al. 1990, 1994). Then 
we select the asteroid with the lowest number (catalogue 
designation) from remaining (not associated) asteroids and 
repeat the HCM association again and again, until no aster- 
oids are left. Then we repeat the whole procedure recursively 
for all clusters detected at Ucutoff = 300 m/s, but now for a 
lower value, e.g., iicutoff = 299 m/s. We may continue until 
Vcutoff = Om/s, but of course, for too low values of the cut- 
off velocity, no clusters can be detected and all asteroids are 
single. The resulting stalactite diagram at Fig. [11] is simply 
the asteroid number (designation) vs Wcutoff plot: a dot at 
a given place is plotted only if the asteroids belongs to a 
cluster of at least max(5, Ai'min(^^cutoft )) bodies. We are not 
interested in clusters with less than 5 members; they are 
most probably random flukes. 

We can see two prominent clusters among Hildas: the 
first one around the asteroid (153) Hilda itself, and the sec- 
ond one around (1911) Schubart. In the remaining part of 
this Section we discuss each of them separately. 

The stalactite diagram constructed in the same way for 
Zhongguos and Griquas is shown in Fig. 1121 No grouping 
seem to be significant enough to be considered an impact- 
generated cluster. This is consistent with the discussion of 
the (ap,ep,sinip) plots in Section [2TT] 



3.1 Schubart family 

The Schubart family can be distinguished from the remain- 
ing population of Hildas on a large range of cutoff velocities: 
from 50 m/s to more than 100 m/s (Fig. Ill|) . It merges with 
the Hilda family at 200 m/s. For the purpose of our analysis 
we selected Ucutoff ~ 60 m/s as the nominal value. While the 
total number of Schubart family members is not too sensitive 
to this cutoff value, we refrain from using too high Vcutoff, for 
which we would expect and increasing number of interlop- 
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Figure 12. A stalactite diagram computed for the long-lived J2/1 
population. There are no prominent groupings; 60 asteroids are 
not associated with any others, even at ticutoff = 300m/s. Every 
group plotted here has at least 5 members. 
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Figure 13. Cumulative distributions N{<H) of absolute mag- 
nitudes H for the whole J3/2 population (solid) and for the 
members of the families around (153) Hilda (dashed) and (1911) 
Schubart (dotted). Thin dotted curves denote N(<H) of the J3/2 
population with one and both families subtracted. Labels are the 
best-fitted values of 7 using N(<H) oc 10''*'^ (straight lines) in the 
interval (12.25, 14.25) for Hilda and (12.25, 14.75) for Schubart 
clusters. For sake of a rough comparison, the upper abscissa gives 
an estimate of the sizes for an albedo value py = 0.044, average 
of the J3/2 population. 



ers to be associated with the family, and the family would 
attain a rather peculiar shape in the (ap, ep, sin ip) space. 

Figure shows the cumulative distribution of the ab- 
solute magnitudes for the Schubart family members, com- 
pared to the rest of the J3/2 population. Importantly, the 
slope 7 = (-f 0.48 ± 0.02) of the N{<H) oc 10^-^ fit is quite 
steeper for the Schubart family, which supports the hypoth- 
esis of its coUisional origin. 



© 2008 RAS, MNRAS OOO.[Tlf20l 



10 M. Broz and D. Vokrouhlicky 




o 

S 15 



Q> 

3 



10 
5 


15 
10 
5 





Hilda 




Schiibart 



1000 



0.1 



0.2 0.3 0.4 0.5 
spectral slope (PC^) 



0.6 0.7 



Figure 14. Distribution of spectral slopes (PCi components of 
the 5 broad-band colours) of 153 asteroids in the J3/2 resonance 
(top), 21 Hilda family members (middle) and 4 Schubart family 
members (bottom). Data computed from the 3rd release of the 
SDSS catalogue of moving objects (ADR3; Ivezic et al. 2002). 
Note some objects have been observed multiple times by SDSS 
and the histograms show distribution of all observations (rather 
than averages for a given body ) . The slope values of Hildas range 
from neutral to steep, with roughly two groups separated by the 
value PCi = 0.3: (i) C/X-types with PCi < 0.3, and (ii) D-types 
with PCi > 0.3. Importantly, the Schubart family members are 
spectrally similar; the median value PCi = 0.20 corresponds to a 
C- or X-type parent body. 



We also analysed the available SDSS catalogue of mov- 
ing objects (ADR3; Ivezic et al. 2002). We searched for the 
J3/2 asteroids among the entries of this catalogue and com- 
puted the principal component PCi of the spectrum in the 
visible band. Note the PCi value is an indicator of the spec- 
tral slope and allows thus to broadly distinguish principal 
spectral classes of asteroids (e.g., Bus et al. 2002). Figure [T^ 
shows our results. The top panel confirms the bimodal char- 
acter of the J3/2 population (see also Dahlgren et al. 1997, 
1998, 1999 and Gil-Hutton & Brunini 2008). More impor- 
tantly, though, the bottom panel indicates a spectral homo- 
geneity of the Schubart family, placing all members within 
the C/X taxonomy class branch. This finding strongly sup- 
ports coUisional origin of the Schubart family. 

Tedesco et al. (2002) derive D = 80 km size for (1911) 
Schubart, corresponding to a very low albedo pv = 0.025. 
The same authors determine D — 38 km size of (4230) van 
den Bergh and exactly the same albedo; this asteroid is 
among the five largest in the family. Assuming the same 
albedo for all other family members, we can construct a size- 
frequency distribution (Fig. I15|) . The slope a ~ (—2. 7 ±0.1) 
fitted to the small end of the distribution, where we still 
assume observational completeness, is rather shallow, but 
marginally within the limits of population slopes produced 
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Figure 15. Cumulative size distribution N{>D) for the whole 
J3/2 population (solid) and for members of the two suggested 
coUisional clusters: Hilda (dashed) and Schubart (dotted). We 
assumed the geometric albedo pv = 0.044 (and 0.025 in the 
Schubart-family case; see the text) for the conversion of absolute 
magnitudes H to diameters D. Labels are the best-fitted values of 
a using N{>D) oc D" (straight lines) in the interval (11, 20) km 
for the two families and (9, 20) km for the J3/2 population. 



in the numerical simulations of disruptions (e.g., Durda et al. 
2007) Q 

If we sum the volumes of the observed members, we end 
up with a lower limit for the parent body size DpB = 110 km, 
provided there are no interlopers. We can also estimate the 
contribution of small (unobserved) bodies using the follow- 
ing simple method: (i) we sum only the volumes of the ob- 
served bodies larger than an assumed completeness limit 
-Dcompicto = 10km (Vcompiotc = Y,ijDf); (ii) we fit the 
cumulative size distribution by a power law (log N{>D) = 
alog[i3]kn, + /3; a = -2.68, /3 = 4.73 for the Schubart); 
(iii) we prolong this slope from Dcompicto down to Dmin = 
and calculate the total volume of the parent body (provided 
a > -3): 



VPB ~ Vcompiotc + — 10 - 

6 a + S 



d: 



a + 3 



complctcj 



(9) 



The result is Dpb = y^Vps = 130 km, some sort of an up- 
per limit. The volumetric ratio between the largest fragment 
and the parent body is then Vlf/Vpb = 0.2, a fairly typical 
value for asteroid families in the main asteroid belt. Obvi- 
ously, the assumption of a single power-law extrapolation 
of the N(>D) at small sizes is only approximate and can 
lead to a result with a 10 % uncertainty. However, if we use 
an entirely different geometric method, developed by Tanga 



^ We also mention that so far asteroid disruption simulations did 
not explore cases of weak-strength materials appropriate for the 
suggested C /X spectral taxonomy of the Schubart family parent 
body. 
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et al. (1998), we obtain Dpb — 120-130 km, i.e., comparable 
to our previous estimate. 

What is an approximate size ddisrupt of a projectile nec- 
essary to disrupt the parent body of the Schubart family? 
Using Eq. (1) from Bottke et al. (2005): 

rfdisrupt = (2Qd/'K^p)^ Dtarget . (10) 

and substituting — 10^ J/kg for the strength (somewhat 
lower than that of basaltic objects to accommodate the as- 
sumed C/X spectral type; e.g., Kenyon et al. 2008 and ref- 
erences therein), Vh-np = 4.78 km/s for the typical impact 
velocity (see Dahlgren 1998) and Dtargot — 130 km, we ob- 
tain ddisrupt — 25 km. At this size the projectile population 
is dominated by main belt bodies. Considering also different 
intrinsic coUisional probabilities between Hilda-Hilda aster- 
oids (2.3 X 10"^* km"^ yr"^; Dahlgren 1998) and Hilda-main 
belt asteroids (0.6 x 10~^* km~^ yr~^), we find it more likely 
the Schubart family parent body was hit by a projectile orig- 
inating from the main belt. 



3.2 Hilda family 

We repeated the same analysis as in Section 13.11 for the 
Hilda family. The family remains statistically distinct from 
the whole J3/2 population in the range of cutoff velocities 
(130, 170) m/s; we choose Wcutoff ~ 150 m/s as the nominal 
value. 

The slope 7 of the cumulative absolute magnitude dis- 
tribution N{<H) is (+0.50 ±0.02) (Fig. [131), again steeper 
than for the total J3/2 population and comparable to that 
of the Schubart family. The spectral slopes (PCi) are some- 
what spread from flat (C/X-compatible values; PCi < 0.3) 
to redder (D-compatible values; PCi > 0.3) — see Fig. 1141 
Overall, though, the C/X members prevail such that the D- 
type objects might be actually interlopers, at least according 
to a simple estimate based on the volume of the Hilda family 
in the (op, Bp, sin ip) space, compared to the total volume of 
the J3/2 population. 

Tedesco et al. (2002) determine albedos for six family 
members. They range from 0.037 to 0.087, but three val- 
ues are close to the median albedo 0.044 of all J3/2 aster- 
oids. We thus consider this value to be representative of the 
Hilda family. The corresponding cumulative size distribution 
is plotted in Fig. 1151 Using the same method as in Section lSTTI 
we estimate the size of the parent body Dpb ~ 180-190 km, 
with Vlf/Vpb ^ 0.8. With the model of Tanga et al. (1999) 
we would obtain Dpb — 210 km and thus Vlf/Vpb — 0.5. 
This family forming event seems to be thus characterized in 
between the catastrophic disruption and a huge cratering. 
The necessary projectile size is ddisrupt ~ 50-55 km. 

While not so prominent as the Schubart family, we con- 
sider the group of asteroids around Hilda a fairly robust case 
of a coUisionally-born family too. 



3.3 Simulated disruption events 

In order to asses some limits for the age of the Schubart and 
Hilda families, we perform a number of numerical tests. In 
particular, we simulate a disruption of a parent body inside 
the resonance and numerically determine the long-term or- 
bital evolution of fragments. The evolved synthetic family 
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Figure 16. The initial osculating elements of an impact- 
generated swarm of 139 fragments at the location of 
(1911) Schubart (bottom crosses), the corresponding pseudo- 
proper elements computed from the first My of evolution (upper 
crosses) and the pseudo-proper elements of the observed Schubart 
family (circles). We show here projection onto the plane defined 
by semimajor axis and eccentricity. Dots are the pseudo-proper 
elements of the background J3/2-population asteroids. The ini- 
tial synthetic swarm of asteroids poorly matches the observed 
family: it is both too extent in semimajor axis and too compact 
in eccentricity. 
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Figure 17. The same as Fig. 1161 but for a diff'erent impact ge- 
ometry (/ = 135°, LU + f = 180°). This choice of / maximizes 
the initial spread of the synthetic family in proper eccentricity. In 
this case we do not show the initial osculating orbital elements. 



at different time steps is then compared with the observed 
family. Ideally, this approach should allow to constrain the 
time elapsed since the family formed. 

As a flrst step, we need to create a synthetic family in- 
side the resonance. We use current orbital elements of the 
largest family member, (1911) Schubart in this case, as rep- 
resentative to the parent body and only allow changes in 
the true anomaly / and in the argument of pericentre oj at 
the break-up event. By changing these two geometric pa- 
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rameters we can produce diflerent initial positions of the 
fragments in the orbital element space. For sake of our test, 
fragments are assumed to be dispersed isotropically with re- 
spect to the parent body, with a velocity distribution given 
by the model of Farinella et al. (1993, 1994). The number 
of fragments dN(v) launched with relative velocities in the 
interval {v, v + dv) is given by 



dN{v) = Cv{v'^ + Wosc 



-{k+1)/2 



dv 



(11) 



with C a normalization constant, Vcsc the escape velocity 
from the parent body and k = 3.25. To prevent excessive 
escape velocities we introduce a maximum allowed value 
"Umax. Nominally, we set Wmax = 200 m/s, but in the next 
Section 13.41 we also use restricted values of this parameter 
to test sensitivity of our results to initial conditions. 

To simulate an impact that might have created the 
Schubart family, we generated velocities randomly for 
139 fragments with v^sc = 65 m/s (note the number of frag- 
ments in the synthetic family is equal to the number of 
the Schubart family members). The resulting swarm of frag- 
ments is shown in Fig. 1161 for the impact geometry / = 0° 
and oj+f = 180°. We show both the initial osculating orbital 
elements and the pseudo-proper elements. 

The synthetic family extends over significantly larger 
range of the semimajor axis than the observed Schubart 
family, but all fragments still fall within the J3/2 resonance. 
The eccentricity distribution is, on the other hand, substan- 
tially more compact. Only the distribution of inclinations of 
the synthetic family roughly matches that of the observed 
family. We verified this holds also for other isotropic-impact 
geometries (such as / = 135° and lu + f = 180° shown in 
Fig. I17|) . The peculiar shape of the synthetic family in the 
pseudo-proper element space {ap,ep) is an outcome of the 
isotropic disruption, simply because some fragments fall to 
the left from the libration centre of the J3/2 resonance (at 
3.97 AU) and they are 'mapped' to the right. This is because 
the pseudo-proper elements are the maxima and minima of 
a and e, respectively, over their resonant oscillations. 

The initial configuration of the synthetic family was 
propagated for 4 Gyr, using the integrator described in Sec- 
tion [2] At this stage, we use only the gravitational perturba- 
tions from the 4 exterior giant planets. We performed such 
simulation for several impact geometries, as determined by 
/ and uj, with similar results. 

Figure [TS] shows the long-term evolution of the syn- 
thetic family. Because the family resides mostly in the sta- 
ble zone of the J3/2 resonance, only little evolution can be 
seen for most of the bodies. This is in accord with findings 
of Nesvorny & Ferraz-Mello (1997) who concluded that the 
stable region in this resonance shows little or no diffusion 
over time scales comparable to the age of the Solar system. 
Only about 10 % of orbits that initially started at the out- 
skirts of the stable zone (with large libration amplitudes) 
escaped from the resonance during the 4 Gyr simulation. 

The removal of orbits with large semimajor axis ap helps 
in part to reconcile the mismatch with the distribution of 
the observed Schubart family. However, the dispersion in 
eccentricity Cp does not evolve much and it still shows large 
mismatch if compared to the observed family. Even in the 
case / = 135° (cj -I- / = 180° ; Fig. I17|l , which maximizes the 
initial eccentricity dispersion of the synthetic fragments, the 
final value at 4 Gyr is about three times smaller than that 
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Figure 18. The synthetic family from Fig. ll6l evolved over 4 Gyr: 
the grey dots show evolutionary tracks of the fragments in the 
pseudo-proper orbital element space. Overall, stability of the J3/2 
resonance makes many fragments to stay very close to their ini- 
tial values. Only ~10% of fragments with the initial extremal 
values of ap (and thus the libration amplitude) escape from the 
resonance during the simulation. This helps in part to reduce the 
mismatch with the observed family (circles) in semimajor axis, 
but is not sufficient to attain the Schubart-family full eccentricity 
dispersion. 



of the Schubart family. Clearly, our model is missing a key 
element to reproduce the current orbital configuration of 
this family. 

One possibility to resolve the problem could be to re- 
lease the assumption of an isotropic impact and explore 
anisotropics in the initial velocity field. This is an obvious 
suspect in all attempts to reconstruct orbital configurations 
of the asteroid families, but we doubt it might help much in 
this case. Exceedingly large relative velocities, compared to 
the escape velocity of the estimated parent body, would be 
required. Recall, the fragments located in the stable region 
of the J3/2 resonance would hardly evolve over the age of 
the Solar system. 

A more radical solution is to complement the force 
model, used for the long-term propagation, by additional 
effects. The only viable mechanism for the size range we are 
dealing with is the Yarkovsky effect. This tiny force, due to 
anisotropic thermal emission, has been proved to have deter- 
mining role in understanding fine structures of the asteroid 
families in the main belt (e.g., Bottke et al. 2001; Vokrouh- 
licky et al. 2006a,b). In these applications the Yarkovsky 
effect produces a steady drift of the semimajor axis, leaving 
other orbital elements basically constant. 

However, the situation is different for resonant orbits. 
The semimajor axis evolution is locked by the strong grav- 
itational influence of Jupiter. For that reason we first ran 
simplified simulations with the Yarkovsky forces — results 
of these tests are briefly described in the Appendix A. We 
next applied the model containing both gravitational and 
Yarkovsky perturbations to the evolution of the synthetic 
family. Results of these experiments are described in the 
next Section [3.41 
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Figure 19. The impact-generated swarm from Fig. 1161 evolved 
by planetary perturbations and the Yarkovsky forces, in the pro- 
jection on the pseudo-proper semimajor axis ap vs eccentricity 
Cp plane. The grey dots indicate the evolutionary tracks over the 
whole 4 Gyr time span and crosses denote the configuration at 
1.7 Gyr, when the eccentricity dispersion of the synthetic family 
particles roughly matches that of the observed Schubart family 
(circles). 
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Figure 20. The difference between the cumulative distribution 
of eccentricities for the observed Schubart family and our syn- 
thetic families expressed as the Kolmogorov-Smirnov (KS) dis- 
tance i?KS vs time t. All models assume / = 0° and u) + f = 180° 
but they have different maximum values t)max of the ejected par- 
ticles: ■Umax = 40, 60, 80, 100 and 200 m/s (labels). For t < 1 Gyr 
and t > 2.4 Gyr the KS distance is large such that the KS prob- 
ability is p(>Dks) < 0.01. Thus with a 99% probability level we 
can exclude such age for the Schubart family. The best matches 
are achieved for Wmax = 40 and 60 m/s, for which the most ex- 
treme particles have been eliminated from the synthetic family. 



3.4 Yarkovsky drift in eccentricity 

We ran our previous simulation of the long-term evolution 
of the synthetic family with the Yarkovsky forces included. 
Our best guess of thermal parameters for bodies of the C/X 
type is: ps = Pb = 1300 kg/m'^ for the surface and bulk densi- 
ties, K = 0.01 W/m/K for the surface thermal conductivity, 
C = 680 J/kg for the heat capacity, A — 0.02 for the Bond 
albedo and e — 0.95 for the thermal emissivity parameter. 
Rotation periods are bound in the 2-12 hours range. Spin 
axes orientations are assumed isotropic in space. Finally, we 
assign sizes to our test particles equal to the estimate of sizes 
for Schubart family members, based on their reported abso- 
lute magnitudes and albedo pv = 0.025. The dependence of 
the Yarkovsky force on these parameters is described, e.g., 
in Bottke et al. (2002, 2006). We note the uncertainties of 
the thermal parameters, assigned to individual bodies, do 
not affect our results significantly, mainly because we simu- 
late a collective evolution of more than 100 bodies; we are 
not interested in evolution of individual orbits. 

We let the synthetic family evolve for 4 Gyr and 
recorded its snapshots every 50 kyr. As discussed in the Ap- 
pendix A, the resonant Yarkovsky effect produces mainly 
secular changes of eccentricity. We recall this systematic 
drift in e must not be confused with the chaotic diffusion 
in e. We also note that inclination of the orbits remains sta- 
ble, in accord with a good match of the Schubart family by 
the initial inclination distribution. 

Because the initial eccentricity dispersion of the syn- 
thetic family is much smaller than that of the observed 
one, its steady increase due to the combined effects of the 
Yarkovsky forces and the resonant lock gives us a possibility 
to date the origin of the family (see Vokrouhlicky et al. 2006a 
for a similar method applied to families in the main belt). 
To proceed in a quantitative way, we use a 1-dimensional 



Kolmogorov-Smirnov (KS) test to compare cumulative dis- 
tribution of pseudo-proper eccentricity values Cp of the ob- 
served and synthetic families (e.g., Press et al. 2007). 

Figure [20] shows the KS distance Dks of the two ec- 
centricity distributions as a function of time. For sake of 
a test, we also use smaller t^max values of the initial veloc- 
ity field (essentially, this is like to start with a more com- 
pact synthetic family). Regardless of the Wmax value, our 
model rejects Schubart family ages smaller than 1 Gyr and 
larger than 2.4 Gyr with a 99 % confidence level. For ages 
in between 1.5 Gyr and 1.7 Gyr the KS-tested likelihood of 
a similarity of the synthetic-family and the observed-family 
Ep distributions can reach up to 50 %. We thus conclude the 
most likely age of the Schubart cluster is (1.7 ± 0.7) Gyr. 

We repeated the same analysis for Hilda family by cre- 
ating a synthetic family of 233 particles about (153) Hilda. 
In this case we used Wosc = 110 m/s. The situation is actu- 
ally very similar to the Schubart — there is again a problem 
with the small dispersion of eccentricities in case of a purely 
gravitational model. Using the model with the Yarkovsky 
effect, we can eventually fit the spread of eccentricities and 
according to the KS test (Fig. I21|l the age of the family 
might be > 4 Gyr. The match is still not perfect, but this 
problem might be partly due to numerous interlopers in the 
family. We also note that a 10 % relative uncertainty of the 
mean albedo of the Hilda family members would lead to a 
5 % uncertainty of their sizes and, consequently, to a 5 % 
uncertainty of the family age. 

The Hilda family seems to be dated back to the Late 
Heavy Bombardment era (Morbidelli et al. 2005; see also 
Sec. 14. ip . We would find such solution satisfactory, because 
the population of putative projectiles was substantially more 
numerous than today (note that a disruption of the Hilda 
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Figure 21. The same as Fig. 1201 but for the Hilda family. The 
KS distance Dks of the observed and modelled pseudo-proper 
eccentricity Cp distribution for several synthetic clusters with dif- 
ferent maximum velocities iimax are plotted vs time t. In this case, 
t ^ 3.5 Gyr seems to be ruled out at a 99 % probability level. 



family parent body is a very unlikely event during the last 
3.5 Gyr). 

We finally simulated a putative collision in the J2/1 res- 
onance, around the asteroid (3789) Zhongguo (the largest 
asteroid in the stable island B). There are two major differ- 
ences as compared to the J3/2 resonance: (i) the underlying 
chaotic diffusion due to the gravitational perturbations is 
larger in the J2/1 resonance (e.g., Nesvorny & Ferraz-Mello 
1997), such that an initially compact cluster would fill the 
whole stable region in 1-1.5 Gyr and consequently becomes 
unobservable; (ii) sizes of the observed asteroids are gener- 
ally smaller, which together with a slightly smaller helio- 
centric distance, accelerates the Yarkovsky drift in e^. The 
latter effect would likely shorten the time scale to 0.5-1 Gyr. 
Thus the non-existence of any significant orbital clusters in 
the J2/1 resonance (Sec. [3] and Fig. I12|) does not exclude 
a coUisional origin of the long-lived population by an event 
older than 1 Gyr. This would also solve the apparent prob- 
lem of the very steep size distribution of the stable J2/1 
population (Broz et al. 2005 and Fig. [3]). Note the expected 
coUisional lifetime of the smallest observed J2/1 asteroids is 
several Gyr (e.g., Bottke et al. 2005). 



4 RESONANT POPULATION STABILITY 
WITH RESPECT TO PLANETARY 
MIGRATION AND THE YARKOVSKY 
EFFECT 

We finally pay a brief attention to the overall stability of as- 
teroid populations in the first-order resonances with respect 
to different configurations of giant planets. We are focus- 
ing on the situations when the orbits of Jupiter and Saturn 
become resonant. This is motivated by currently adopted 
views about final stages of building planetary orbits archi- 
tecture, namely planet migration in a diluted planetesimal 
disk (e.g., Malhotra 1995; Hahn & Malhotra 1999; Tsiga- 
nis et al. 2005; MorbideUi et al. 2005; Gomes et al. 2005 - 
these last three references are usually described as the Nice 
model). MorbideUi et al. (2005) proved that the primordial 
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Figure 22. Dependence of the period associated with circulation 
of the Great Inequality (GI) angle 5Ag — 2Aj (Ag and Aj are 
mean longitude in orbit for Saturn and Jupiter) on the osculating 
semimajor axis ag of Saturn (crosses); Jupiter's orbit is fixed. 
Moving Saturn farther away from the location of 2:5 resonance 
with Jupiter (at ag = 9.584 AU) makes the GI period shorter 
(with today's ag = 9.5289 AU the GI period is about 880 yr). 
In the interval of ag values shown in this figure the GI period 
becomes comparable to the average period Pj2/i of librations in 
the J2/1 resonance (solid line). Libration period for a particular 
case of asteroid (1362) Griqua is shown by circles. 



Trojan asteroids were destabilized when Jupiter and Sat- 
urn crossed their mutual 1:2 mean motion resonance and, 
at the same time, the Trojan region was re-populated by 
particles of the planetesimal disk. Since the mutual 1:2 res- 
onance of Jupiter and Saturn plays a central role in the 
Nice model, and since these two planets had to cross other 
(weaker) mutual resonances such as 4:9 and 3:7 before they 
acquired today's orbits, one can naturally pose a question 
about the stability of primordial populations in the first- 
order mean motion resonances with Jupiter. Ferraz-Mello 
et al. (1998a,b) demonstrated that even subtler effects can 
influence the J2/1 population, namely the resonances be- 
tween the asteroid libration period in the J2/1 resonance and 
the period of the Great Inequality (GI) terms in planetary 
perturbations (i.e., those associated with Jupiter and Sat- 
urn proximity to their mutual 2:5 mean motion resonance; 
Fig. I22|l . A first glimpse to the stability of the first-order 
resonance populations with respect to these effects is given 
in Section [4.11 

In Section 14.21 we also briefly estimate the change of 
dynamical lifetimes for small J2/1 and J3/2 bodies caused 
by the Yarkovsky effect. 



4.1 Planetary migration effects 

In what follows we use a simple approach by only moving 
Saturn's orbit into different resonance conflgurations with 
Jupiter's orbit. We do not let orbits of these planets migrate, 
but consider them static. With such a crude approach we can 
only get a flrst hint about a relative role of depletion of the 
asteroid populations in the flrst-order resonances (note in 
reality the planets undergo steady, but likely not smooth, 
migration and exhibit jumps over different mutual resonant 
states; e.g., supplementary materials of Tsiganis et al. 2005). 
The results are summarised in Fig. 1231 (i) The Hilda 
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Figure 23. Histograms of dynamical lifetimes for asteroids in the 
J2/1 [top], J3/2 [middle] and J4/3 [bottom] resonances, in case 
Jupiter and Saturn are at their current orbits, or in a mutual 1:2, 
4:9, 3:7, 2:5 resonance, or in a Great Inequality resonance (in case 
of the J2/1 only). The histograms were computed for 106 long- 
lived asteroids in the J2/1, first 100 Hildas in the J3/2 and 8 in 
the J4/3 (including short-lived). 



group in the J3/2 resonance is very unstable (on the time 
scale ~ 1 Myr) with respect to the 1:2 Jupiter-Saturn reso- 
nancejfl on the contrary J2/1 asteroids may survive several 
10 Myr in this configuration of Jupiter and Saturn, so this 
population is not much affected by this phase by planetary 
evolution (note, moreover, that Jupiter and Saturn would 
likely cross the zone of other mutual 1:2 resonance in ~ 1 Myr 
only; e.g., Tsiganis et al. 2005 and MorbideUi et al. 2005). 
(ii) The 4:9 resonance has a larger infiuence on the J2/1 
population than on Hildas, (iii) In the case of 3:7 resonance 
it is the opposite: the J3/2 is more unstable than the J2/1. 
(iv) The Great Inequality resonance does indeed destabilise 
the J2/1 on a time scale 50 Myr. Provided the last phases of 



° Jupiter Trojans, which are already known to be strongly unsta- 
ble (MorbideUi et al. 2005), would have the dynamical lifetime of 
the order 0.1 Myr in this kind of simulation. 
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Figure 24. An A'^-body simulation of planetary migration driven 
by dynamically cold planetesimal disk beyond Neptune, with 10^ 
particles and total mass 50 , and including also 10^ mass-less 
particles in the J3/2 resonance with Jupiter. Top: The semima- 
jor axis as of Saturn vs time and the position of the 1:2 mean 
motion resonance with Jupiter (estimated from the Kepler law 
(]^/2)-2/3aj). Bottom: The same for asteroids initially located 
inside the J3/2 resonance with Jupiter. The J3/2 asteroids are 
strongly destabilised at the very time of the 1:2 Jupiter-Saturn 
resonance crossing [t = 1.25 Myr) and none of the 1000 test parti- 
cles survived in the J3/2 region after a mere 0.5 Myr. This means 
more than 99.9 % depletion of the primordial population. None of 
the planetesimal-disk particles got trapped in the J3/2 during or 
after Jupiter's and Saturn's passage through the 1:2 resonance, 
indicating that more particles are needed to study this process. 
We used the Mercury hybrid-scheme integrator (Chambers 1999) 
for the purpose of this test. The gravitational interactions be- 
tween planets and massive planetesimals are accounted for, but 
planetesimals do not interact with each other, nor with mass-less 
test particles. The time step was 36 days and the accuracy pa- 
rameter 10"^'^. Initial conditions of planets were: aj = 5.2 AU, 
ag = 8.05 AU, axj = 12.3 AU, = 17.5 AU, with all eccentrici- 
ties and inclination of the order 10~^. We took the current orbits 
of Hildas as the initial conditions for our test particles. Note the 
destabilisation of the Hilda region is neither sensitive to precise 
initial conditions nor to the mass of the planetesimal disk; the 
only relevant condition is that Jupiter and Saturn cross their mu- 
tual 1:2 resonance. 



the migration proceed very slowly, it may cause a significant 
depletion of the primordial J2/1 population. In the exact 2:5 
resonance, the J2/1 population would not be affected at all. 

According to our preliminary tests with a more com- 
plete Ai'-body model for planetary migration which includes 
a disk of 10'' planetesimals beyond Neptune, the strong 
instability of the J3/2 asteroids indeed occurs during the 
Jupiter-Saturn 1:2 resonance crossing (see Fig. I24|l . A vast 
reservoir of planetesimals residing beyond the giant plan- 
ets, and to some extent also nearby regions of the outer 
asteroid belt, are probably capable to re-populate the J3/2 
resonance zone at the same time and form the currently ob- 
served populations. This is similar to the Trojan clouds of 
Jupiter (MorbideUi et al. 2005). 
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Figure 25. Logarithmic histograms of dynamical lifetimes for 
the originally long-lived asteroids in the J2/1 resonance, in case 
the Yarkovsky effect perturbs the orbits. The sizes of the objects 
(4-12 km) were multiplied by 2, 0.2, 0.02 and 0.002 for compari- 
son. Note a stronger instability starts to occur for the factor 0.02 
(i.e., for the sizes 0.08-0.24 km). The YORP reorientations are 
not included in this model. 



4.2 The Yarkovsky effect 



In Section 13.41 we already discussed the influence of the 
Yarkovsky effect on the families located inside the J3/2 res- 
onance. In course of time it modified eccentricities of their 
members, but did not cause a large-scale instability; the fam- 
ilies remained inside the resonance aU the time. Here we 
seek the size threshold for which the Yarkovsky would case 
overall instability by quickly removing the bodies from the 
resonance. 

We perform the following numerical test: we multiply 
sizes of the long-lived J2/1 objects by fudge factors 0.2, 
0.02 and 0.002, for which the Yarkovsky effect is stronger, 
and compare respective dynamical lifetimes with the original 
long-lived objects. Results are summarized in Fig. 1251 We 
can conclude that a significant destabilisation of the J2/1 
resonant population occurs for sizes ~ 0.1 km and smaller 
(provided the nominal population has sizes mostly 4-12 km). 

We do not include the YORP effect (i.e., the torque in- 
duced by the infrared thermal emission) at this stage. The 
YORP is nevertheless theoretically capable to significantly 
decelerate (or accelerate) the rotation rate, especially of the 
smallest asteroids, which can lead to random reorientations 
of the spin axes due to collisions, because the angular mo- 
mentum is low in this spin-down state. These reorientations 
can be simulated by a Monte-Carlo model with a typical 
time scale (Capek & Vokrouhlicky 2004; R is radius in kilo- 
metres and a orbital semimajor axis in AU): 



TYORP ^ 25 Myr {R/l kmf (2.5 A\J/af 



(12) 



Since the Yarkovsky effect depends on the obliquity value, 
the systematic drift would be changed to a random walk for 
bodies whose spin axis undergo frequent re-orientations by 
the YORP effect. We can thus expect that the YORP effect 
might significantly prolong dynamical lifetimes of resonant 
objects with sizes ~ 0.1 km or smaller, because tyorp < 
0.25 Myr for them. 

We can also check, which orientation of the spin axis 
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Figure 26. The same as Fig. 1251 but now the Yarkovsky effect is 
calculated only for small bodies (0.08— 0.24 km) and varies due to 
different values of the obliquity 7 = 0, 45°, 90°, 135° and 180°. 
Note the retrograde bodies (i.e., with a negative drift da/dt < 0) 
are significantly more unstable. 



makes the escape from the J2/1 resonance more likely 
to happen. We consider 0.08-0.24 km bodies, clone them 
5 times and assign them different values of the obliquity 
-y = 0°, 45°, 90°, 135° and 180°. Figure shows clearly 
that the retrograde rotation increases the probability of the 
escape. This is consistent with the structure of the J2/1 
resonance, for which low-a separatrix does not continue to 
e = Q. We conclude the remaining very small (yet unobserv- 
able) asteroids inside the J2/1 may exhibit a preferential 
prograde rotation. 

We perform a similar simulation for the J3/2 population 
(first 100 bodies with sizes 10-60 km), but now with YORP 
reorientations (Eq. I12p included. The results (Fig. I27p show 
the J3/2 population is much less affected than the J2/1 by 
the Yarkovsky/ YORP perturbation. 

We conclude the Yarkovsky/YORP effect may destabi- 
lize the J2/1 and J3/2 bodies only partially on the 100 Myr 
time scale and only for sizes smaller than ~ 0.1 km. It is ob- 
viously a remote goal to verify this conclusion by observa- 
tions (note the smallest bodies in these resonances have sev- 
eral kilometers size). Nevertheless, the dynamical lifetimes 
of small asteroids determined in this section, might be useful 
for collisional models of asteroid populations, which include 
also dynamical removal. 



5 CONCLUSIONS AND FUTURE WORK 

The main results of this paper can be summarised as fol- 
lows: (i) we provided an update of the observed J2/1, J3/2 
and J4/3 resonant populations; (ii) we discovered two new 
objects in the J4/3 resonance; (iii) we described two as- 
teroid families located inside the J3/2 resonance (Schubart 
and Hilda) and provided an evidence that they are of a col- 
lisional origin; (iv) we reported a new mechanism how the 
Yarkovsky effect systematically changes eccentricities of res- 
onant asteroids; we used this phenomenon to estimate the 
ages of the Schubart and Hilda families ((1.7 ± 0.7) Gyr and 
> 4 Gyr respectively); (v) coUisionally-born asteroid clus- 
ters in the stable region of J2/1 would disperse in about 
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Figure 27. The same as Fig. 1251 but for the J3/2 resonance. 
InstabiUty occurs for the sizes multiplied by 0.002 (i.e., 0.02- 
0.12 km). The YORP effect and corresponding reorientations of 
the spin axes are included in this case. 
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Figure 28. Pseudo-proper eccentricity Cp vs absolute magnitude 
H plot for the members of the Hilda family. The triangular shape 
and the depletion of asteroids in the centre of the family (around 
Bp ~ 0.22) is discussed in the text. 



1 Gyr; (vi) 20% of Hildas may escape from the J3/2 reso- 
nance within 4 Gyr in the current configuration of planets; 
(vii) Hildas are strongly unstable when Jupiter and Saturn 
cross their mutual 1:2 mean motion resonance. 

The J3/2 resonance is a unique 'laboratory' — the 
chaotic diffusion is so weak, that families almost do not dis- 
perse in eccentricity and inclination due to this effect over 
the age of the Solar system. What is even more important, 
they almost do not disperse in semimajor axis, even thought 
the Yarkovsky effect operates. The drift in a is transformed 
to a drift in e, due to a strong gravitational coupling with 
Jupiter. We emphasize, this is is not a chaotic diffusion in e, 
but a systematic drift in e. 

Another piece of information about the families in J3/2 
resonance is hidden in the eccentricity Bp vs absolute magni- 
tude H plots (see Fig.[28]for the Hilda family). The triangu- 
lar shape (larger eccentricity dispersion of the family mem- 
bers for larger H) is a well-known combination of two ef- 
fects: (i) larger ejection speed, and (ii) faster dispersal by the 
Yarkovky forces for smaller fragments. Interestingly, there 
is also a noticeable depletion of small bodies in the centre 
of the family and their concentration at the outskirts — a 
phenomenon known from (a, H) plots of main-belt families, 
which was interpreted as an interplay between the Yarkovsky 
and YORP effects (Vokrouhhcky et al. 2006)0 Indeed the 
estimated ~ 4 Gyr age for this family matches the time scale 
of a YORP cycle for D ~ 10 km asteroids in the Hilda region 
(e.g., Vokrouhlicky & Capek 2002; Capek & Vokrouhlicky 
2004) . Vokrouhlicky et al. (2006a) pointed out that this cir- 
cumstance makes the uneven distribution of family members 
most pronounced. 

We postpone the following topics for the future work: 
(i) a more precise age determination for the resonant aster- 
oid families, based on the Yarkovsky /YORP evolution in the 
(e, H) space; (ii) a more detailed modelling of analytic or A''- 



^ The YORP effect tilts the spin axes of asteroids preferentially 
towards obliquity 7 = 0° or 180° and this enhances the diurnal 
Yarkovsky drift due to its cos 7 dependence. 



body migration of planets and its influence on the stability 
of resonant populations. 



ACKNOWLEDGEMENTS 

We thank David Nesvorny, Alessandro Morbidelli and 
William F. Bottke for valuable discussions on the subject 
and also Rodney Gomes for a constructive review. The work 
has been supported by the Grant Agency of the Czech Re- 
public (grants 205/08/0064 and 205/08/P196) and the Re- 
search Program MSM0021620860 of the Czech Ministry of 
Education. We also acknowledge the usage of the Metacen- 
trum computing facilities (http://meta.cesnet.cz/) and 
computers of the Observatory and Planetarium in Hradec 
Kralove (http : //www . astrohk .cz/[ . 



REFERENCES 

Beauge C, Roig F., 2001, Icarus, 153, 391 

Bottke W. F., Vokrouhlicky D., Broz M., Nesvorny D., Morbidelli 

A., 2001, Sci, 294, 1693 
Bottke W. F., Vokrouhlicky D., Rubincam D. P., Broz M., 2002, 

in Bottke W.F., Cellino A., Paolicchi P., Binzel R.P., eds.. 

Asteroids III. University of Arizona Press, Tucson, p. 395 
Bottke W. F., Durda D. D., Nesvorny D., Jedicke R., Morbidelli 

A., Vokrouhlicky D., Levison H. F., 2005, Icarus, 179, 63 
Bottke W. F., Vokrouhlicky D., Rubincam D. P., Nesvorny D., 

2006, Annu. Rev. Earth Planet. Sci., 34, 157 
Brown M.E., Barkume K.M., Ragozzinc D., Schaller E.L., 2007, 

Nature, 446, 294 
Broz M., Vokrouhlicky D., Roig F., Nesvorny D., Bottke W.F., 

Morbidelli A., 2005, MNRAS, 359, 1437 
Broz M., 2006, PhD thesis, Charles University (see at 

http:/ /sirrah. troja.mff.cuni.cz/~mira/mp/) 
Bus S.J., Vilas F., Barucci M.A., 2002, in Bottke W.F., CeUino 

A., Paolicchi P., Binzel R.P., eds.. Asteroids III. University of 

Arizona Press, Tucson, p. 169 
Capek D., Vokrouhlicky D., 2004, Icarus, 172, 526 
Chambers I.E., 1999, MNRAS, 304, 793 
Dahlgren M., 1998, A&A, 336, 1056 
Dahlgren M., Lagerkvist C.-I., 1995, A&A, 302, 907 



© 2008 RAS, MNRAS OOO.[Tlf20l 



18 M. Broz and D. Vokrouhlicky 



Dahlgren M., Lagerkvist C.-L, Fitzsimmons A., Williams I. P., 
Gordon M., 1997, A&A, 323, 606 

Dahlgren M., Lahulla J.F., Lagerkvist C.-I., Lagerros J., Mot- 
tola S., Erikson A., Gonano-Beurer M., di Martino M., 1998, 
Icarus, 133, 247 

Dahlgren M., Lahulla J.F., Lagerkvist C.-L, 1999, Icarus, 138, 259 
Davis D.R., Durda D.D., Marzari F., Campo Bagatin A., Gil- 
Hutton R., 2002, in Bottke W.F., CeUino A., PaoUcchi P., 
Binzel R.P., eds., Asteroids III. University of Arizona Press, 
Tucson, p. 545 

DeirOro A., Marzari F., PaoUcchi P., Vanzani A., 2001, A&A, 
366, 1053 

Dohnanyi J.W., 1969, J. Geophys. Res., 74, 2531 

Durda D., Bottke W.F., Nesvorny D., Encke B., Merline W.J., 

Asphaug E., Richardson D.C., 2007, Icarus, 186, 498 
Farinella P., Gonczi R., Froeschle Ch., Froeschle C, 1993, Icarus, 

101, 174 

Farinella P., Froeschle C, Gonczi R., 1994, in Milani A., Di Mar- 
tino M., CeUino A., eds.. Asteroids, comets, meteors 1993. 
Kluwer Academic Publishers, Dordrecht, p. 205 

Ferraz-Mello S., 1988, AJ, 96, 400 

Ferraz-Mello S., Michtchenko T.A., 1996, A&A, 310, 1021 
Ferraz-Mello S., Michtchenko T.A., Nesvorny D., Roig F., Simula, 

A., 1998a, Planet. Sp. Sci., 46, 1425 
Ferraz-Mello S., Michtchenko T.A., Roig F., 1998b, AJ, 116, 1491 
Gil-Hutton R., Brunini A., 2008, Icarus, 193, 567 
Gomes R., Levison H. L., Tsiganis K., Morbidelli A., 2005, Na- 
ture, 435, 466 
Hahn J. M., Malhotra R., 1999, AJ, 117, 3041 
Henrard J., 1982, Celest. Mech., 27, 3 

Ivezic Z., Juric M., Lupton R.H., Tabachnik S., Quinn T., 2002, 
Survey and other telescope technologies and discoveries, (J. A. 
Tyson, S. Wolff, Eds.) Proceedings of SPIE Vol. 4836 

Jedicke R., Magnier E.A., Kaiser N., Chambers K.C., 2007, in Mi- 
lani A., Valsecchi G.B., Vokrouhlicky D., eds.. Near Earth As- 
teroids, our Celestial Neighbors. Cambridge University Press, 
Cambridge, p. 341 

Kenyon S. J., Bromley B. C, O'Brien D. P., Davis D. R., 2008, in 
Barucci M.A., Boehnhardt H., Cruikshank D.P., Morbidelli 
A., eds.. The Solar System beyond Neptune. University of 
Arizona Press, Tucson, p. 293 

Krueger A., 1889, AN, 120, 221 

Kuhnert F., 1876, AN, 87, 173 

Landau L. D., Lifschiz E. M., 1976, Mechanics, Pergamon Press, 
Oxford 

Laskar J., Robutel P., 2001, Celest. Mech. Dyn. Astron., 80, 39 
Lemaitre A., Henrard J., 1990, Icarus, 83, 391 
Levison H., Duncan M., 1994, Icarus, 108, 18 
Malhotra R., 1995, AJ, 110, 420 

Michtchenko T. A., Ferraz-Mello S., 1997, Planet. Sp. Sci., 45, 
1587 

Milani A., 1993, Celest. Mech. Dyn. Astron., 57, 59 
Milani A., FarineUa P., 1995, Icarus, 115, 209 
MorbidelU A., Moons M., 1993, Icarus, 102, 316 
MorbidelU A., Levison H. L., Tsiganis K., Gomes R., 2005, Na- 
ture, 435, 462 

Morbidelli A., Tsiganis K., Crida A., Levison H.F., Gomes R., 

2007, AJ, 134, 1790 
Moons M., MorbidelU A., Migliorini F., 1998, Icarus, 135, 458 
Murray C. D., 1986, Icarus, 65, 70 

Murray C. D., Dermott S. F., 1999, Solar System Dynamics, Cam- 
bridge University Press, Cambridge 
Nesvorny D., Ferraz-Mello S., 1997, Icarus, 130, 247 
Nesvorny D., Beauge C, Dones L., 2004, AJ, 127, 1768 
Nesvorny D., Vokrouhlicky D., Bottke W.F., 2006, Science, 312, 
1490 

Nesvorny D., Bottke W.F., Dones L., Levison H.F., 2002, Nature, 
417, 720 



Nesvorny D., Alvarellos J.L.A., Dones L., Levison H.F., 2003, AJ, 
126, 398 

Nesvorny D., Jedicke R., Whiteley R.J., Ivezic Z., 2005, Icarus, 
173, 132 

O'Brien D.P., Greenberg R., 2003, Icarus, 164, 334 
Palisa J., 1888, AN, 120, 111 

Press W. R., Teukolsky S. A., Vetterling W., Flannery B. P., 

2007, Numerical Recipes: The Art of Scientific Computing, 

Cambridge University Press, Cambridge 
Quinn T.R., Tremaine, S., Duncan, M. 1991, AJ, 101, 2287 
Roig F., Gil-Hutton R., 2006, Icarus, 183, 411 
Roig F., Nesvorny D., Ferraz-Mello S., 2002, MNRAS, 335, 417 
Roig F., Ribeiro A. O., Gil-Hutton R., 2008, A&A, in press 
Schubart J., 1982a, A&A, 114, 200 
Schubart J., 1982b, Cel. Mech. Dyn. Astron., 28, 189 
Schubart J., 1991, A&A, 241, 297 
Schubart J., 2007, Icarus, 188, 189 
Sessin W., Bressane B.R., 1988, Icarus, 75, 97 
Tanga P., CeUino A., Michel P., Zappala V., PaoUcchi P., Dell'Oro 

A., 1998, Icarus, 141, 65 
Tedesco E.F., Noah P.V., Noah M., Price S.D., 2002, AJ, 123, 

1056 

Tsiganis K., Knezevic Z., VarvogUs H., 2007, Icarus, 186, 484 
Tsiganis K., Gomes R., MorbidelU A., Levison H.L., 2005, Nature, 
435, 459 

Tsuchida M., 1990, Rev. Mexicana Astron. Astrophys., 21, 585 

Vokrouhlicky D., Broz, 2002, in Celletti A., Ferraz-Mello S., Hen- 
rard J., eds.. Modern Celestial Mechanics: from Theory to 
Applications. Kluwer Academic Publishers, Dordrecht, p. 467 

Vokrouhlicky D., Broz M., Bottke W.F., Nesvorny D., Morbidelli 
A., 2006a, Icarus, 182, 118 

Vokrouhlicky D., Broz M., MorbidelU A., Bottke W.F., Nesvorny 
D., Lazzaro D., Rivkin A. S., 2006b, Icarus, 182, 92 

Weissman P.R., Bottke W.F., Levison H.L., 2002, in Bottke W.F., 
CeUino A., PaoUcchi P., Binzel R.P., eds.. Asteroids HI. Uni- 
versity of Arizona Press, Tucson, p. 669 

Zappala V., CeUino A., Farinella P., Knezevic Z., 1990, AJ, 100, 
2030 

Zappala V., CeUino A., Farinella P., Milani A., 1994, AJ, 107, 772 
Zappala V., CeUino A., DeU'Oro A., PaoUcchi P., 2002, in Bottke 

W.F., CeUino A., PaoUcchi P., Binzel R.P., eds.. Asteroids HI. 

University of Arizona Press, Tucson, p. 619 



APPENDIX A: RESONANT YARKOVSKY 
EFFECT 

The effects of weak dissipative forces, such as the tidal force, 
gas-drag force and the Poynting-Robertson force, on both 
non-resonant and resonant orbits were extensively studied 
in the past (e.g., Murray & Dermott 1999 and references 
therein). Interaction of the Yarkovsky drifting orbits with 
high-order, weak resonances was also numerically studied to 
some extent (e.g., Vokrouhlicky & Broz 2002) but no sys- 
tematic effort was paid to study Yarkovsky evolving orbits 
in strong low-order resonances. Here we do not intend to 
develop a detailed theory, rather give a numerical example 
that can both help to explain results presented in the main 
text and motivate a more thorough analytical theory. 

The Yarkovsky effect outside the resonance. We con- 
structed the following simple numerical experiment: we took 
the current orbit of (1911) Schubart as a starting condition 
and integrated the motion of two 0.1km size objects with 
extreme obliquity values 0° and 180° . Their thermal param- 
eters were the same as in Section 13.41 Because the diurnal 
variant of the Yarkovsky effect dominates the evolution, the 
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Figure Al. Orbital evolution of two D = 0.1 km asteroids in the J3/2 resonance within a circular restricted three-body problem (Jupiter 
on a circular orbit). Obliquities 0° and 180° were assigned to the two bodies, such that outside the resonance they would migrate by the 
Yarkovsky forces in opposite direction. The expected change Aa of the semimajor axis in 100 My is depicted by the arrow on top of the 
left panel. The left and middle panels show 1 kyr orbital segments at the beginning and at the end of the simulation: (i) in the semimajor 
axis a vs eccentricity e projection (left), and (ii) in the projection of Cartesian resonant variables \/2S (cos a, sin a) (see Eqs. [T] and (2] 
short-period variations have been removed for better visibility) (middle). The orbits slowly evolve from the initial A^o — 0.45 level-curve 
of the integral given by Eq. ([5J to their final values of ~ 0.44 (Af+ with da/dt > 0) and ~ 0.46 (A'^_ with da/dt < 0), respectively (see 
also the right panel). During this evolution the libration centre follows the position of the exact periodic orbit in the J3/2 (dashed curve 
in the left panel). Because the latter has a steep progression in e as a changes, orbital evolution is characterised by a significant change 
of the eccentricity Ae (also Bp) but only a small change in a (also ap). 



extreme obliquities would mean the two test bodies would 
normally (outside any resonances) drift in semimajor axis 
in two opposite directions (e.g., Bottke et al. 2002, 2006). 
The two orbits would secularly acquire Aa ~ -1-0.25 AU or 
—0.25 AU in 100 My, about the extent shown by the arrow 
on top of the left panel of Fig. lAll Since the strength of the 
Yarkovsky forces is inversely proportional to the size, we can 
readily scale the results for larger bodies. 

The resonance without the Yarkovsky effect. If we in- 
clude gravitational perturbations by Jupiter only, within a 
restricted circular three body problem (e,j = 0), and remove 
short-period oscillations by a digital filter, the parameter 
A'^ from Eq. (O would stay constant. The orbit would be 
characterized by a stable libration in (S,cr) variables with 
about 30° amplitude in a (see the curve labelled Myr in 
the middle panel of Fig. IA1|) . 

While evolving, some parameters known as the adia- 
batic invariants are approximately conserved (see, e.g.. Lan- 
dau & Lifschitz 1976; Henrard 1982; Murray & Dermott 
1999). One of the adiabatic invariants is N itself. Another, 
slightly more involved quantity, is the area J enclosed by 
the resonant path in the V2E (cos ct, sin a) space: 

J = ^V2Eda. (Al) 

We would thus expect these parameters be constant, except 
for strong-enough perturbation or long-enough time scales 
(recall the adiabatic invariants are constant to the second 
order of the perturbing parameter only). 

Resonant Yarkovsky effect. Introducing the Yarkovsky 
forces makes the system to evolve slowly. The lock in the res- 
onance prevents the orbits to steadily drift away in the semi- 
major axis and the perturbation by the Yarkovsky forces 
acts adiabatically. This is because (i) the time scale of the 
resonance oscillations is much shorter than the characteristic 
time scale of the orbital evolution driven by the Yarkovsky 



forces, and (ii) the strength of the resonant terms in the 
equations of motion are superior to the strength of the 
Yarkovsky accelerations. 

We let the two J3/2 orbits evolve over 100 Myr 
(Fig. lAlfl . At the end of our simulation the orbits moved 
from A'^o ~ 0.45 to N+ ~ 0.46 (for the outward migrating 
orbit) or to N- ~ 0.44 (for the inward migrating orbit), re- 
spectively. During this evolution, both orbits remained per- 
manently locked in the J2/1 resonance, librating about the 
periodic orbit (dashed line in the left panel of Fig. lAip . Be- 
cause the position of this centre has a steep progression in 
the eccentricity and only small progression in the semima- 
jor axis, the evolution across different A-planes makes the 
orbital eccentricity evolve significantly more than the semi- 
major axis. This is also seen in the middle panel of Fig. lAll 
where the librating orbits significantly split farther/closer 
with respect to origin of coordinates (note the polar dis- 
tance from the origin is basically a measure of the eccen- 
tricity). The shape of the librating orbit is modified such 
that the area J stays approximately constant. We have ver- 
ified that the relative change in both adiabatic invariants, 
acquired during the 100 Myr of evolution, is about the same: 
5N/N ~ SJ/J ~ 5 X lO"'^. It is a direct expression of the 
strength of the perturbation by the Yarkovsky forces. 

We can conclude the Yarkovsky eff'ect results in a sig- 
nificantly different type of secular evolution for orbits ini- 
tially inside strong first-order mean motion resonances with 
Jupiter. Instead of secularly pushing the orbital semimajor 
axis inward or outward from the Sun, it drives the orbital 
eccentricity to smaller or larger values, while leaving the 
semimajor axis to follow the resonance centre. 

If we were to leave the orbital evolution continue in 
our simple model, the inward-migrating orbit would leave 
the resonance toward the zone of low-eccentricity apocen- 
tric librators. Such bodies are observed just below the J2/1 
resonance. On the other hand, the outward-migrating orbit 
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would finally increase the eccentricity to the value when the 
orbit starts to cross the Jupiters orbit. Obviously, in a more 
complete model, with all planets included, the orbits would 
first encounter the unstable region surrounding the stable 
resonant zone. Such marginally stable populations exist in 
both the J3/2 and J2/1 resonances. 

This paper has been typeset from a TJ^X/ I^Tj^X file prepared 
by the author. 
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